
# Replication code for: 
# "Who Benefits? Perceived Policy Effects and Beneficiaries of Agricultural Protection in Japan"
# Hirofumi Kawaguchi (Keio SFC) & Ikuma Ogura (Hitotsubashi)
# Last updated on January 15, 2026

# Preparation----

# Load packages
require(dplyr)
require(estimatr)
require(ggplot2)
require(haven)
require(marginaleffects)
require(MASS)
require(mediation)
require(modelsummary)
options("modelsummary_format_numeric_latex" = "plain")

# Load data 
dat <- read_dta("JJES_replication.dta")

# Create directory to store figures
if (sum(grepl("Figure", list.files())) == 0) dir.create("./Figure")

# Variable labels for regression tables
var.labs <- c("(Intercept)" = "(Intercept)",
              treatDevelopment = "Development of Agriculture",
              treatDisaster = "Disaster Prevention",
              "treatFood sufficiency" = "Food Sufficiency",
              treatQuality = "Quality of Agriculture",
              treatRural = "Richer Rural Communities")

# Seed value for random number generation
set.seed(1234)

# Session information
sessionInfo()

# Data preprocessing----

# Experimental variables
dat$FMC.true <- dplyr::case_when(
  grepl("食料自給率", dat$merit) ~ 1,
  grepl("持続的", dat$merit) ~ 2,
  grepl("豊か", dat$merit) ~ 3,
  grepl("品質", dat$merit) ~ 4,
  grepl("災害", dat$merit) ~ 5,
  TRUE ~ 8
)
dat$FMC.correct <- (dat$FMC == dat$FMC.true) * 1 # FMC
prop.table(table(dat$FMC.correct))

dat$treat <- dplyr::case_when(
  grepl("食料自給率", dat$merit) ~ "Food sufficiency",
  grepl("持続的", dat$merit) ~ "Development",
  grepl("豊か", dat$merit) ~ "Rural",
  grepl("品質", dat$merit) ~ "Quality",
  grepl("災害", dat$merit) ~ "Disaster",
  TRUE ~ "Control"
)
dat$target <- dplyr::case_when(
  grepl("周りの人", dat$for_whom) ~ "Extended",
  grepl("日本社会", dat$for_whom) ~ "Sociotropic",
  TRUE ~ "Pocketbook"
)
dat$dv1 <- 6 - as.numeric(dat$Experiment_1)
dat$dv2 <- ifelse(dat$Experiment_2 == 1, 3, 
                  ifelse(dat$Experiment_2 == 2, 1, 2))

# Other variables
dat$age_cat <- factor(dat$Age,
                      levels = c(1:5), 
                      labels = c("20-29", "30-39", "40-49", "50-59", "60-69"))
dat$gender <- factor(dat$Sex,
                     levels = c(1, 2, 3),
                     labels = c("男", "女", "その他"))
dat$region <- factor(dat$Prefecture_1,
                     levels = c(1, 3, 10, 18, 28, 36, 42, 47, 56),
                     labels = c("北海道", "東北", "関東", "中部", 
                                "近畿", "中国", "四国", "九州", "その他"))
dat$city <- ifelse(dat$Citysize == 1, "Large", 
                   ifelse(dat$Citysize == 2, "Medium",
                          ifelse(dat$Citysize == 3, "Small", NA))) |> factor()
dat$college <- (dat$Educ == 5 | dat$Educ == 6) * 1
dat$DQ.correct <- ifelse(dat$Issue_3 == 1 & dat$Issue_8 == 4, 1, 0) 
dat$primary.industry <- (dat$Industry == 1) * 1
dat$agri.know <- ifelse(dat$Know_1 == 1 | dat$Know_2 == 1 | dat$Know_3 == 1, 1, 0)

# Main analysis: Analysis 1----
out1.1 <- lm_robust(dv1 ~ treat 
                    + age_cat + gender + city
                    + college + factor(Income) + as.numeric(Ideology_1)
                    + agri.know + primary.industry, 
                    data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                         & target == "Sociotropic"))
out1.2 <- lm_robust(dv1 ~ treat 
                    + age_cat + gender + city
                    + college + factor(Income) + as.numeric(Ideology_1)
                    + agri.know + primary.industry, 
                    data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                         & target == "Extended"))
out1.3 <- lm_robust(dv1 ~ treat 
                    + age_cat + gender + city
                    + college + factor(Income) + as.numeric(Ideology_1)
                    + agri.know + primary.industry, 
                    data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                         & target == "Pocketbook"))
modelsummary(list("Sociotropic" = out1.1, 
                  "Extended Self-interest" = out1.2, 
                  "Self-interest" = out1.3), 
             stars = TRUE, output = "latex",
             coef_map = var.labs,
             gof_omit = "AIC|BIC|RMSE") # Table A1
# Plot it
tiff("./Figure/Fig2.tiff", 
     width = 6.5, height = 4, units = "in", compression = "lzw", res = 1200)
layout(matrix(c(1, 2, 3), ncol = 3), widths = c(2.5, 1, 1))
par(mar = c(3, 0.5, 3, 0.5))
## Sociotropic
plot(NULL, NULL, type = "n", 
     xlim = c(-1.25, 0.5), ylim = c(0.5, 5.5),
     xlab = "", ylab = "", axes = FALSE)
text(-1.2, c(5:1), 
     labels = c("Development of Agriculture", "Disaster Prevention", "Food Sufficiency",
                "Quality Agricultural Products", "Richer Rural Communities"),
     pos = 4, cex = 1.2)
text(0, 5.7, substitute(paste(bold("Sociotropic"))), xpd = NA, cex = 1.2)
tmp.col <- ifelse(confint(out1.1)[2:6, 1] > 0, 1, "gray50")
segments(0, 0.5, 0, 5.5, col = "gray80")
segments(-0.15, c(5:1), 0.45, c(5:1), col = "gray80", lty = 2)
points(coef(out1.1)[2:6], c(5:1), pch = 19, col = tmp.col)
segments(confint(out1.1)[2:6, 1], c(5:1), confint(out1.1)[2:6, 2], c(5:1),
         col = tmp.col)
axis(1, at = seq(-0.15, 0.45, 0.15))
## Extended self-interest
plot(NULL, NULL, type = "n", 
     xlim = c(-0.2, 0.5), ylim = c(0.5, 5.5),
     xlab = "", ylab = "", axes = FALSE)
text(0, 5.7, substitute(paste(bold("Extended Self-Interest"))), xpd = NA, cex = 1.2)
tmp.col <- ifelse(confint(out1.2)[2:6, 1] > 0, 1, "gray50")
segments(0, 0.5, 0, 5.5, col = "gray80")
segments(-0.15, c(5:1), 0.45, c(5:1), col = "gray80", lty = 2)
points(coef(out1.2)[2:6], c(5:1), pch = 19, col = tmp.col)
segments(confint(out1.2)[2:6, 1], c(5:1), confint(out1.2)[2:6, 2], c(5:1),
         col = tmp.col)
axis(1, at = seq(-0.15, 0.45, 0.15))
## Extended self-interest
plot(NULL, NULL, type = "n", 
     xlim = c(-0.2, 0.5), ylim = c(0.5, 5.5),
     xlab = "", ylab = "", axes = FALSE)
text(0, 5.7, substitute(paste(bold("Self-Interest"))), xpd = NA, cex = 1.2)
tmp.col <- ifelse(confint(out1.3)[2:6, 1] > 0, 1, "gray50")
segments(0, 0.5, 0, 5.5, col = "gray80")
segments(-0.15, c(5:1), 0.45, c(5:1), col = "gray80", lty = 2)
points(coef(out1.3)[2:6], c(5:1), pch = 19, col = tmp.col)
segments(confint(out1.3)[2:6, 1], c(5:1), confint(out1.3)[2:6, 2], c(5:1),
         col = tmp.col)
axis(1, at = seq(-0.15, 0.45, 0.15))
dev.off()

# Main analysis: Analysis 2----
dat$Treat <- (dat$treat != "Control") * 1
# Sociotropic
m.1 <- lm(dv1 ~ Treat + age_cat + gender + city
          + college + factor(Income) + as.numeric(Ideology_1)
          + agri.know, 
          data = dat |> subset(DQ.correct == 1 & Industry != 1
                               & target == "Sociotropic"))
y.1 <- lm(dv2 ~ Treat + dv1 + age_cat + gender + city
          + college + factor(Income) + as.numeric(Ideology_1)
          + agri.know, 
          data = dat |> subset(DQ.correct == 1 & Industry != 1
                               & target == "Sociotropic"))
my.1 <- mediate(m.1, y.1, treat = "Treat", mediator = "dv1")
summary(my.1)
# Extended self-interest
m.2 <- lm(dv1 ~ Treat + age_cat + gender + city
          + college + factor(Income) + as.numeric(Ideology_1)
          + agri.know, 
          data = dat |> subset(DQ.correct == 1 & Industry != 1
                               & target == "Extended"))
y.2 <- lm(dv2 ~ Treat + dv1 + age_cat + gender + city
          + college + factor(Income) + as.numeric(Ideology_1)
          + agri.know, 
          data = dat |> subset(DQ.correct == 1 & Industry != 1
                               & target == "Extended"))
my.2 <- mediate(m.2, y.2, treat = "Treat", mediator = "dv1")
summary(my.2)
# Self-interest
m.3 <- lm(dv1 ~ Treat + age_cat + gender + city
          + college + factor(Income) + as.numeric(Ideology_1)
          + agri.know, 
          data = dat |> subset(DQ.correct == 1 & Industry != 1
                               & target == "Pocketbook"))
y.3 <- lm(dv2 ~ Treat + dv1 + age_cat + gender + city
          + college + factor(Income) + as.numeric(Ideology_1)
          + agri.know, 
          data = dat |> subset(DQ.correct == 1 & Industry != 1
                               & target == "Pocketbook"))
my.3 <- mediate(m.3, y.3, treat = "Treat", mediator = "dv1")
summary(my.3)
# Plot it
tiff("./Figure/Fig3.tiff", 
     width = 6.5, height = 4, units = "in", compression = "lzw", res = 1200)
par(mar = c(3, 0.5, 3, 0.5))
plot(NULL, NULL, type = "n", 
     xlim = c(-0.2, 0.12), ylim = c(0.5, 3.5),
     axes = FALSE, xlab = "", ylab = "")
text(-0.2, c(3:1), 
     labels = c("Sociotropic Evaluation", "Extended Self-interest Evaluation", "Self-interest Evaluation"),
     pos = 4, cex = 1.2)
text(-0.16, 3.7, substitute(paste(bold("Average Causal Mediation Effects via..."))), 
     xpd = NA, pos = 4, cex = 1.2)
segments(0, 0.5, 0, 3.5, col = "gray80")
segments(-0.02, c(3:1), 0.12, c(3:1), col = "gray80", lty = 2)
points(c(my.1$d1, my.2$d1, my.3$d1), c(3:1), pch = 19, col = c(1, "gray50", 1))
segments(c(my.1$d1.ci[1], my.2$d1.ci[1], my.3$d1.ci[1]), c(3:1), 
         c(my.1$d1.ci[2], my.2$d1.ci[2], my.3$d1.ci[2]), c(3:1), 
         col = c(1, "gray50", 1))
axis(1, at = seq(-0.02, 0.12, 0.02))
dev.off()

# Main analysis: Analysis 3----
# Experimental treatment -> Outcome 1
citysize <- c("Large", "Medium", "Small")
out.list1 <- vector("list", length = 9)
for (i in 1:3){
  dat.sub <- dat |> subset(DQ.correct == 1 & Industry != 1 & city == citysize[i]) 
  out.list1[[(i - 1) * 3 + 1]] <- lm_robust(dv1 ~ treat 
                                            + age_cat + gender
                                            + college + factor(Income) + as.numeric(Ideology_1) 
                                            + agri.know, 
                                            data = dat.sub |> subset(target == "Sociotropic"))
  out.list1[[(i - 1) * 3 + 2]] <- lm_robust(dv1 ~ treat 
                                            + age_cat + gender
                                            + college + factor(Income) + as.numeric(Ideology_1) 
                                            + agri.know, 
                                            data = dat.sub |> subset(target == "Extended"))
  out.list1[[(i - 1) * 3 + 3]] <- lm_robust(dv1 ~ treat 
                                            + age_cat + gender
                                            + college + factor(Income) + as.numeric(Ideology_1) 
                                            + agri.know, 
                                            data = dat.sub |> subset(target == "Pocketbook"))
}
modelsummary(list("Sociotropic" = out.list1[[1]],
                  "Extended Self-interest" = out.list1[[2]],
                  "Self-interest" = out.list1[[3]]),
             stars = TRUE, output = "latex",
             coef_map = var.labs,
             gof_omit = "AIC|BIC|RMSE")
modelsummary(list("Sociotropic" = out.list1[[4]],
                  "Extended Self-interest" = out.list1[[5]],
                  "Self-interest" = out.list1[[6]]),
             stars = TRUE, output = "latex",
             coef_map = var.labs,
             gof_omit = "AIC|BIC|RMSE")
modelsummary(list("Sociotropic" = out.list1[[7]],
                  "Extended Self-interest" = out.list1[[8]],
                  "Self-interest" = out.list1[[9]]),
             stars = TRUE, output = "latex",
             coef_map = var.labs,
             gof_omit = "AIC|BIC|RMSE") # Table A3
tiff("./Figure/Fig4.tiff",
     width = 6.5, height = 9, units = "in", compression = "lzw", res = 1200)
layout(matrix(c(1:9), nrow = 3, ncol = 3, byrow = TRUE), widths = c(2.5, 1, 1))
par(mar = c(3, 0.5, 3, 0.5))
for (i in 1:3){
  # Sociotropic
  plot(NULL, NULL, type = "n", 
       xlim = c(-2.43, 0.9), ylim = c(0.5, 5.5),
       xlab = "", ylab = "", axes = FALSE)
  text(-2.4, c(5:1), 
       labels = c("Development of Agriculture", "Disaster Prevention", "Food Sufficiency",
                  "Quality Agricultural Products", "Richer Rural Communities"),
       pos = 4, cex = 1.2)
  text(0, 5.7, substitute(paste(bold("Sociotropic"))), xpd = NA, cex = 1.2)
  segments(0, 0.5, 0, 5.5, col = "gray80")
  segments(-0.45, c(5:1), 0.9, c(5:1), col = "gray80", lty = 2)
  points(coef(out.list1[[(i - 1) * 3 + 1]])[2:6], c(5:1), pch = 16, 
         col = ifelse(confint(out.list1[[(i - 1) * 3 + 1]])[2:6, 1] > 0, 1, "gray50"))
  segments(confint(out.list1[[(i - 1) * 3 + 1]])[2:6, 1], c(5:1), confint(out.list1[[(i - 1) * 3 + 1]])[2:6, 2], c(5:1),
           col = ifelse(confint(out.list1[[(i - 1) * 3 + 1]])[2:6, 1] > 0, 1, "gray50"))
  axis(1, at = seq(-0.4, 0.8, 0.2), cex.axis = 0.9)
  if (i == 1){
    legend(-2.7, 6, substitute(paste(bold("(a) Large City"))), bty = "n", xpd = TRUE, cex = 1.2)
  } else if (i == 2) {
    legend(-2.7, 6, substitute(paste(bold("(b) Mid-sized City"))), bty = "n", xpd = TRUE, cex = 1.2)
  } else {
    legend(-2.7, 6, substitute(paste(bold("(c) Small City/Town/Village"))), bty = "n", xpd = TRUE, cex = 1.2)
  }
  # Extended self-interest
  plot(NULL, NULL, type = "n", 
       xlim = c(-0.45, 0.9), ylim = c(0.5, 5.5),
       xlab = "", ylab = "", axes = FALSE)
  text(0, 5.7, substitute(paste(bold("Extended Self-Interest"))), xpd = NA, cex = 1.2)
  segments(0, 0.5, 0, 5.5, col = "gray80")
  segments(-0.45, c(5:1), 0.9, c(5:1), col = "gray80", lty = 2)
  points(coef(out.list1[[(i - 1) * 3 + 2]])[2:6], c(5:1), pch = 16, 
         col = ifelse(confint(out.list1[[(i - 1) * 3 + 2]])[2:6, 1] > 0, 1, "gray50"))
  segments(confint(out.list1[[(i - 1) * 3 + 2]])[2:6, 1], c(5:1), confint(out.list1[[(i - 1) * 3 + 2]])[2:6, 2], c(5:1),
           col = ifelse(confint(out.list1[[(i - 1) * 3 + 2]])[2:6, 1] > 0, 1, "gray50"))
  axis(1, at = seq(-0.4, 0.8, 0.2), cex.axis = 0.9)
  # Self-interest
  plot(NULL, NULL, type = "n", 
       xlim = c(-0.45, 0.9), ylim = c(0.5, 5.5),
       xlab = "", ylab = "", axes = FALSE)
  text(0, 5.7, substitute(paste(bold("Self-Interest"))), xpd = NA, cex = 1.2)
  segments(0, 0.5, 0, 5.5, col = "gray80")
  segments(-0.45, c(5:1), 0.9, c(5:1), col = "gray80", lty = 2)
  points(coef(out.list1[[(i - 1) * 3 + 3]])[2:6], c(5:1), pch = 16, 
         col = ifelse(confint(out.list1[[(i - 1) * 3 + 3]])[2:6, 1] > 0, 1, "gray50"))
  segments(confint(out.list1[[(i - 1) * 3 + 3]])[2:6, 1], c(5:1), confint(out.list1[[(i - 1) * 3 + 3]])[2:6, 2], c(5:1),
           col = ifelse(confint(out.list1[[(i - 1) * 3 + 3]])[2:6, 1] > 0, 1, "gray50"))
  axis(1, at = seq(-0.4, 0.8, 0.2), cex.axis = 0.9)
}
dev.off()
# Mediation analysis
out.mat1 <- matrix(NA, nrow = 9, ncol = 3)
out.list.med1 <- vector("list", length = 9)
for (i in 1:3){
  dat.sub <- dat |> subset(DQ.correct == 1 & primary.industry == 0 & city == citysize[i]) 
  m.1 <- lm(dv1 ~ Treat + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1)
            + agri.know, 
            data = dat.sub |> subset(target == "Sociotropic"))
  y.1 <- lm(dv2 ~ Treat + dv1 + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1)
            + agri.know, 
            data = dat.sub |> subset(target == "Sociotropic"))
  out.list.med1[[(i - 1) * 3 + 1]] <- tmp.1 <- mediate(m.1, y.1, treat = "Treat", mediator = "dv1")
  m.2 <- lm(dv1 ~ Treat + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1)
            + agri.know, 
            data = dat.sub |> subset(target == "Extended"))
  y.2 <- lm(dv2 ~ Treat + dv1 + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1)
            + agri.know, 
            data = dat.sub |> subset(target == "Extended"))
  out.list.med1[[(i - 1) * 3 + 2]] <- tmp.2 <- mediate(m.2, y.2, treat = "Treat", mediator = "dv1")
  m.3 <- lm(dv1 ~ Treat + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1)
            + agri.know, 
            data = dat.sub |> subset(target == "Pocketbook"))
  y.3 <- lm(dv2 ~ Treat + dv1 + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1)
            + agri.know, 
            data = dat.sub |> subset(target == "Pocketbook"))
  out.list.med1[[(i - 1) * 3 + 3]] <- tmp.3 <- mediate(m.3, y.3, treat = "Treat", mediator = "dv1")
  out.mat1[(i - 1) * 3 + 1,] <- c(tmp.1$d1, tmp.1$d1.ci)
  out.mat1[(i - 1) * 3 + 2,] <- c(tmp.2$d1, tmp.2$d1.ci)
  out.mat1[(i - 1) * 3 + 3,] <- c(tmp.3$d1, tmp.3$d1.ci)
}
summary(out.list.med1[[1]]); summary(out.list.med1[[2]]); summary(out.list.med1[[3]])
summary(out.list.med1[[4]]); summary(out.list.med1[[5]]); summary(out.list.med1[[6]])
summary(out.list.med1[[7]]); summary(out.list.med1[[8]]); summary(out.list.med1[[9]]) # Table A4
tiff("./Figure/Fig5.tiff",
     width = 6.5, height = 4, units = "in", compression = "lzw", res = 1200)
layout(matrix(c(1, 2, 3), ncol = 3), widths = c(2.86, 1, 1))
par(mar = c(3, 0.5, 3, 0.5))
# Large City
plot(NULL, NULL, type = "n", 
     xlim = c(-0.8, 0.2), ylim = c(0.5, 3.5),
     xlab = "", ylab = "", axes = FALSE)
text(-0.8, c(3:1), 
     labels = c("Sociotropic Evaluation", "Extended Self-interest Evaluation", "Self-interest Evaluation"),
     pos = 4, cex = 1.1)
text(-0.8, 3.3, substitute(paste(bold("Average Causal Mediation Effect via..."))),
     pos = 4, cex = 1.1)
text(0, 3.7, substitute(paste(bold("Large City"))), xpd = NA, cex = 1.1)
tmp.col <- ifelse(out.mat1[1:3, 2] >= 0 | out.mat1[1:3, 3] <= 0, 1, "gray50")
segments(0, 0.5, 0, 3.5, col = "gray80")
segments(-0.15, c(3:1), 0.2, c(3:1), col = "gray80", lty = 2)
points(out.mat1[1:3, 1], c(3:1), pch = 16, col = tmp.col)
segments(out.mat1[1:3, 2], c(3:1), out.mat1[1:3, 3], c(3:1),
         col = tmp.col)
axis(1, at = seq(-0.1, 0.2, 0.1))
# Mid-sized city
plot(NULL, NULL, type = "n", 
     xlim = c(-0.15, 0.2), ylim = c(0.5, 3.5),
     xlab = "", ylab = "", axes = FALSE)
text(0, 3.7, substitute(paste(bold("Mid-sized City"))), xpd = NA, cex = 1.1)
tmp.col <- ifelse(out.mat1[4:6, 2] >= 0 | out.mat1[4:6, 3] <= 0, 1, "gray50")
segments(0, 0.5, 0, 3.5, col = "gray80")
segments(-0.15, c(3:1), 0.2, c(3:1), col = "gray80", lty = 2)
points(out.mat1[4:6, 1], c(3:1), pch = 16, col = tmp.col)
segments(out.mat1[4:6, 2], c(3:1), out.mat1[4:6, 3], c(3:1),
         col = tmp.col)
axis(1, at = seq(-0.1, 0.2, 0.1))
# Small city
plot(NULL, NULL, type = "n", 
     xlim = c(-0.15, 0.2), ylim = c(0.5, 3.5),
     xlab = "", ylab = "", axes = FALSE)
text(0, 3.7, substitute(paste(bold("Small City/Town/Village"))), xpd = NA, cex = 1.1)
tmp.col <- ifelse(out.mat1[7:9, 2] >= 0 | out.mat1[7:9, 3] <= 0, 1, "gray50")
segments(0, 0.5, 0, 3.5, col = "gray80")
segments(-0.15, c(3:1), 0.2, c(3:1), col = "gray80", lty = 2)
points(out.mat1[7:9, 1], c(3:1), pch = 16, col = tmp.col)
segments(out.mat1[7:9, 2], c(3:1), out.mat1[7:9, 3], c(3:1),
         col = tmp.col)
axis(1, at = seq(-0.1, 0.2, 0.1))
dev.off()

# Additional analysis: Distribution of the outcome variables by city size (Online Appendix C.1)----
dist.out1 <- dist.out2 <- vector("list", length = 3)
for (i in 1:3){
  dat.sub <- subset(dat, treat == "Control" & city == citysize[i])
  dist.out1[[i]] <- tapply(dat.sub$dv1, dat.sub$target, function(x) prop.table(table(x)))
  dist.out2[[i]] <- tapply(dat.sub$dv2, dat.sub$target, function(x) prop.table(table(x)))
}
# Plot it
dv.labs <- c("Sociotropic", "Extended Self-interest", "Self-interest")
cairo_pdf("./Figure/FigA1.pdf", width = 6.5, height = 5, pointsize = 9,
          family = "Helvetica")
par(mfrow = c(3, 3), mar = c(3, 4, 3, 1))
for (i in 1:3){
  for (j in 1:3){
    barplot(dist.out1[[i]][[j]], ylim = c(0, 0.5),
            main = ifelse(i == 1, dv.labs[j], ""), cex.main = 1.5)
    if (j == 1 & i == 1){
      mtext(substitute(paste(bold("Large City"))), side = 2, line = 2)
    } else if (j == 1 & i == 2) {
      mtext(substitute(paste(bold("Mid-sized City"))), side = 2, line = 2)
    } else if (j == 1 & i == 3) {
      mtext(substitute(paste(bold("Small City/Town/Village"))), side = 2, line = 2)
    } 
  }
}
dev.off()
cairo_pdf("./Figure/FigA2.pdf", width = 6.5, height = 5, pointsize = 9,
          family = "Helvetica")
par(mfrow = c(3, 3), mar = c(3, 4, 3, 1))
for (i in 1:3){
  for (j in 1:3){
    barplot(dist.out2[[i]][[j]], ylim = c(0, 0.5),
            main = ifelse(i == 1, dv.labs[j], ""), cex.main = 1.5)
    if (j == 1 & i == 1){
      mtext(substitute(paste(bold("Large City"))), side = 2, line = 2)
    } else if (j == 1 & i == 2) {
      mtext(substitute(paste(bold("Mid-sized City"))), side = 2, line = 2)
    } else if (j == 1 & i == 3) {
      mtext(substitute(paste(bold("Small City/Town/Village"))), side = 2, line = 2)
    } 
  }
}
dev.off()

# Additional analysis: Sensitivity analysis of the mediation analysis (Online Appendix C.2)----
med1 <- medsens(my.1, rho.by = 0.1, effect.type = "indirect", sims = 100)
med2 <- medsens(my.2, rho.by = 0.1, effect.type = "indirect", sims = 100)
med3 <- medsens(my.3, rho.by = 0.1, effect.type = "indirect", sims = 100)

cairo_pdf("./Figure/FigA3.pdf", width = 6.5, height = 4, pointsize = 9,
          family = "Helvetica")
par(mfrow = c(1, 3), mar = c(4, 4, 6, 1))
plot(med1, main = "Sociotropic Evaluation", cex.main = 1.4, cex.lab = 1.4)
plot(med2, main = "Extended Self-interest Evaluation", cex.main = 1.4, cex.lab = 1.4)
plot(med3, main = "Self-interest Evaluation", cex.main = 1.4, cex.lab = 1.4)
mtext(substitute(paste(bold("Average Causal Mediation Effect via..."))), side = 3, outer = TRUE, 
      line = -2, adj = 0, cex = 1.2)
dev.off()

# Additional analysis: Subsample analysis by acquaintance with farmers (Online Appendix C.3)----
# Experimental treatment -> Outcome 1
out.list2 <- vector("list", length = 6)
for (i in 1:2){
  dat.sub <- dat |> subset(DQ.correct == 1 & Industry != 1 & agri.know == (i - 1)) 
  out.list2[[(i - 1) * 3 + 1]] <- lm(dv1 ~ treat 
                                     + age_cat + gender
                                     + college + factor(Income) + as.numeric(Ideology_1), 
                                     data = dat.sub |> subset(target == "Sociotropic"))
  out.list2[[(i - 1) * 3 + 2]] <- lm(dv1 ~ treat 
                                     + age_cat + gender
                                     + college + factor(Income) + as.numeric(Ideology_1), 
                                     data = dat.sub |> subset(target == "Extended"))
  out.list2[[(i - 1) * 3 + 3]] <- lm(dv1 ~ treat 
                                     + age_cat + gender
                                     + college + factor(Income) + as.numeric(Ideology_1), 
                                     data = dat.sub |> subset(target == "Pocketbook"))
}
modelsummary(list("Sociotropic" = out.list2[[1]],
                  "Extended Self-interest" = out.list2[[2]],
                  "Self-interest" = out.list2[[3]]),
             stars = TRUE, output = "latex",
             vcov = "HC1",
             coef_map = var.labs,
             gof_omit = "AIC|BIC|RMSE|Log.Lik.")
modelsummary(list("Sociotropic" = out.list2[[4]],
                  "Extended Self-interest" = out.list2[[5]],
                  "Self-interest" = out.list2[[6]]),
             stars = TRUE, output = "latex",
             vcov = "HC1",
             coef_map = var.labs,
             gof_omit = "AIC|BIC|RMSE|Log.Lik.") # Table A5
# Mediation analysis
out.mat2 <- matrix(NA, nrow = 6, ncol = 3)
out.list.med2 <- vector("list", length = 6)
for (i in 1:2){
  dat.sub <- dat |> subset(DQ.correct == 1 & primary.industry == 0 & agri.know == (i - 1)) 
  m.1 <- lm(dv1 ~ Treat + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1), 
            data = dat.sub |> subset(target == "Sociotropic"))
  y.1 <- lm(dv2 ~ Treat + dv1 + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1), 
            data = dat.sub |> subset(target == "Sociotropic"))
  out.list.med2[[(i - 1) * 3 + 1]] <- tmp.1 <- mediate(m.1, y.1, treat = "Treat", mediator = "dv1")
  m.2 <- lm(dv1 ~ Treat + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1), 
            data = dat.sub |> subset(target == "Extended"))
  y.2 <- lm(dv2 ~ Treat + dv1 + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1), 
            data = dat.sub |> subset(target == "Extended"))
  out.list.med2[[(i - 1) * 3 + 2]] <- tmp.2 <- mediate(m.2, y.2, treat = "Treat", mediator = "dv1")
  m.3 <- lm(dv1 ~ Treat + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1), 
            data = dat.sub |> subset(target == "Pocketbook"))
  y.3 <- lm(dv2 ~ Treat + dv1 + age_cat + gender
            + college + factor(Income) + as.numeric(Ideology_1), 
            data = dat.sub |> subset(target == "Pocketbook"))
  out.list.med2[[(i - 1) * 3 + 3]] <- tmp.3 <- mediate(m.3, y.3, treat = "Treat", mediator = "dv1")
  out.mat2[(i - 1) * 3 + 1,] <- c(tmp.1$d1, tmp.1$d1.ci)
  out.mat2[(i - 1) * 3 + 2,] <- c(tmp.2$d1, tmp.2$d1.ci)
  out.mat2[(i - 1) * 3 + 3,] <- c(tmp.3$d1, tmp.3$d1.ci)
}
summary(out.list.med2[[1]]); summary(out.list.med2[[2]]); summary(out.list.med2[[3]])
summary(out.list.med2[[4]]); summary(out.list.med2[[5]]); summary(out.list.med2[[6]]) # Table A6
# Treatment effect heterogeneity by city size among respondents w/ farmer acquaintances
out.list3 <- vector("list", length = 3)
dat.sub <- dat |> subset(DQ.correct == 1 & Industry != 1 & agri.know == 0) 
out.list3[[1]] <- lm_robust(dv1 ~ Treat*city 
                            + college + factor(Income) + as.numeric(Ideology_1), 
                            data = dat.sub |> subset(target == "Sociotropic"))
out.list3[[2]] <- lm_robust(dv1 ~ Treat*city 
                            + age_cat + gender
                            + college + factor(Income) + as.numeric(Ideology_1), 
                            data = dat.sub |> subset(target == "Extended"))
out.list3[[3]] <- lm_robust(dv1 ~ Treat*city 
                            + age_cat + gender
                            + college + factor(Income) + as.numeric(Ideology_1), 
                            data = dat.sub |> subset(target == "Pocketbook"))
modelsummary(list("Sociotropic" = out.list3[[1]],
                  "Extended Self-interest" = out.list3[[2]],
                  "Self-interest" = out.list3[[3]]),
             stars = TRUE, output = "latex",
             coef_map = c("(Intercept)" = "(Intercept)",
                          Treat = "Experimental Treatment",
                          cityMedium = "Mid-sized City",
                          citySmall = "Small City/Town/Village",
                          "Treat:cityMedium" = "Treatment x Mid-sized City",
                          "Treat:citySmall" = "Treatment x Small City/Town/Village"),
             gof_omit = "AIC|BIC|RMSE") # Table A7
## compute the marginal effects + plot
out.mat3 <- rbind.data.frame(avg_slopes(out.list3[[1]], variables = "Treat", by = "city") |> 
                               subset(select = c("city", "estimate", "conf.low", "conf.high")),
                             avg_slopes(out.list3[[2]], variables = "Treat", by = "city") |> 
                               subset(select = c("city", "estimate", "conf.low", "conf.high")),
                             avg_slopes(out.list3[[3]], variables = "Treat", by = "city") |> 
                               subset(select = c("city", "estimate", "conf.low", "conf.high")))
cairo_pdf("./Figure/FigA4.pdf", width = 6.5, height = 4, pointsize = 9,
          family = "Helvetica")
layout(matrix(c(1, 2, 3), nrow = 1, ncol = 3), widths = c(2, 1, 1))
par(mar = c(3, 1, 3, 1))
plot(NULL, NULL, type = "n", xlim = c(-1.15, 0.55), ylim = c(0.5, 3.5),
     xlab = "", ylab = "", axes = FALSE)
text(0, 3.7, substitute(paste(bold("Sociotropic"))), xpd = NA, cex = 1.5)
text(-1.15, 3:1, c("Large City", "Mid-sized City", "Small City/Town/Village"),
     pos = 4, cex = 1.5)
abline(v = 0, col = "gray80")
segments(-0.3, 3:1, 0.55, 3:1, lty = 2, col = "gray80")
points(out.mat3$estimate[1:3], 3:1, pch = 16,
       col = ifelse(out.mat3$conf.low[1:3] >= 0 | out.mat3$conf.high[1:3] <= 0, 1, "gray50"))
segments(out.mat3$conf.low[1:3], 3:1, out.mat3$conf.high[1:3], 3:1,
         col = ifelse(out.mat3$conf.low[1:3] >= 0 | out.mat3$conf.high[1:3] <= 0, 1, "gray50"))
axis(1, at = seq(-0.2, 0.4, 0.2))
plot(NULL, NULL, type = "n", xlim = c(-0.3, 0.55), ylim = c(0.5, 3.5),
     xlab = "", ylab = "", axes = FALSE)
text(0, 3.7, substitute(paste(bold("Extended Self-interest"))), xpd = NA, cex = 1.5)
abline(v = 0, col = "gray80")
segments(-0.3, 3:1, 0.55, 3:1, lty = 2, col = "gray80")
points(out.mat3$estimate[4:6], 3:1, pch = 16,
       col = ifelse(out.mat3$conf.low[4:6] >= 0 | out.mat3$conf.high[4:6] <= 0, 1, "gray50"))
segments(out.mat3$conf.low[4:6], 3:1, out.mat3$conf.high[4:6], 3:1,
         col = ifelse(out.mat3$conf.low[4:6] >= 0 | out.mat3$conf.high[4:6] <= 0, 1, "gray50"))
axis(1, at = seq(-0.2, 0.4, 0.2))
plot(NULL, NULL, type = "n", xlim = c(-0.3, 0.55), ylim = c(0.5, 3.5),
     xlab = "", ylab = "", axes = FALSE)
text(0, 3.7, substitute(paste(bold("Self-interest"))), xpd = NA, cex = 1.5)
abline(v = 0, col = "gray80")
segments(-0.3, 3:1, 0.55, 3:1, lty = 2, col = "gray80")
points(out.mat3$estimate[7:9], 3:1, pch = 16,
       col = ifelse(out.mat3$conf.low[7:9] >= 0 | out.mat3$conf.high[7:9] <= 0, 1, "gray50"))
segments(out.mat3$conf.low[7:9], 3:1, out.mat3$conf.high[7:9], 3:1,
         col = ifelse(out.mat3$conf.low[7:9] >= 0 | out.mat3$conf.high[7:9] <= 0, 1, "gray50"))
axis(1, at = seq(-0.2, 0.4, 0.2))
dev.off()

# Additional analysis: Other pre-registered analysis (Online Appendix C.4)----
# Experimental treatment -> Outcome 2
out2.1 <- lm_robust(dv2 ~ treat 
                    + age_cat + gender + city
                    + college + factor(Income) + as.numeric(Ideology_1)
                    + agri.know + primary.industry, 
                    data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                         & target == "Sociotropic"))
out2.2 <- lm_robust(dv2 ~ treat 
                    + age_cat + gender + city
                    + college + factor(Income) + as.numeric(Ideology_1)
                    + agri.know + primary.industry, 
                    data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                         & target == "Extended"))
out2.3 <- lm_robust(dv2 ~ treat 
                    + age_cat + gender + city
                    + college + factor(Income) + as.numeric(Ideology_1)
                    + agri.know + primary.industry, 
                    data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                         & target == "Pocketbook"))
modelsummary(list("Sociotropic" = out2.1,
                  "Extended Self-interest" = out2.2,
                  "Self-interest" = out2.3),
             stars = TRUE, output = "latex",
             coef_map = var.labs,
             gof_omit = "AIC|BIC|RMSE") # Table A8

# Ordered logistic regression analyses
## Experimental treatments -> Outcome 1
out3.1 <- polr(ordered(dv1) ~ treat 
               + age_cat + gender + city
               + college + factor(Income) + as.numeric(Ideology_1)
               + agri.know + primary.industry, 
               data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                    & target == "Sociotropic"),
               Hess = TRUE, method = "logistic")
out3.2 <- polr(ordered(dv1) ~ treat 
               + age_cat + gender + city
               + college + factor(Income) + as.numeric(Ideology_1)
               + agri.know + primary.industry, 
               data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                    & target == "Extended"),
               Hess = TRUE, method = "logistic")
out3.3 <- polr(ordered(dv1) ~ treat 
               + age_cat + gender + city
               + college + factor(Income) + as.numeric(Ideology_1)
               + agri.know + primary.industry, 
               data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                    & target == "Pocketbook"),
               Hess = TRUE, method = "logistic")
modelsummary(list("Sociotropic" = out3.1,
                  "Extended Self-interest" = out3.2,
                  "Self-interest" = out3.3),
             stars = TRUE, output = "latex",
             coef_map = var.labs,
             gof_omit = "R2|R2 adj.|RMSE|Log.Lik.|BIC") # Table A9
## Experimental treatment -> Outcome 2
out4.1 <- polr(ordered(dv2) ~ treat 
               + age_cat + gender + city
               + college + factor(Income) + as.numeric(Ideology_1)
               + agri.know + primary.industry, 
               data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                    & target == "Sociotropic"),
               Hess = TRUE, method = "logistic")
out4.2 <- polr(ordered(dv2) ~ treat 
               + age_cat + gender + city
               + college + factor(Income) + as.numeric(Ideology_1)
               + agri.know + primary.industry, 
               data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                    & target == "Extended"),
               Hess = TRUE, method = "logistic")
out4.3 <- polr(ordered(dv2) ~ treat 
               + age_cat + gender + city
               + college + factor(Income) + as.numeric(Ideology_1)
               + agri.know + primary.industry, 
               data = dat |> subset(DQ.correct == 1 & Industry != 1 
                                    & target == "Pocketbook"),
               Hess = TRUE, method = "logistic")
modelsummary(list("Sociotropic" = out4.1,
                  "Extended Self-interest" = out4.2,
                  "Self-interest" = out4.3),
             stars = TRUE, output = "latex",
             coef_map = var.labs,
             gof_omit = "R2|R2 adj.|RMSE|Log.Lik.|BIC") # Table A10

